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I.  SUMMARY 


GENERALIZED  NON-LINEAR  MINIMAL  RESIDUAL  (GNLMR) 

METHOD  FOR  OPTIMAL  MULTISTEPPING 

George  S.  Dulikravich  Penn  State  University 

Department  of  Aerospace  Engineering 
233  Hammond  Building 
University  Park,  PA  16802 

One  of  the  extrapolation  methods  for  the  acceleration  of  iterative  algorithms  is  the 
Generalized  Non-Linear  Minimal  Residual  (GNLMR)  concept.  It  utilizes  a  number  of 
intermediate  steps  when  advancing  the  solution  to  the  next  time  level.  That  is, 
numerical  error  at  the  new  time  level  can  be  expressed  as  a  sum  of  intermediate 
corrections  where  each  correction  is  multiplied  by  a  separate  acceleration  factor  which 
GNLMR  optimizes.  The  method  was  originally  developed  by  Kennon  and  Dulikravich  and 
then  successfully  generalized  and  applied  by  Huang  and  Dulikravich  to  a  number  of 
problems  governed  by  single  non-linear  partial  differential  equations.  In  addition,  Huang 
has  obtained  preliminary  results  where  GNLMR  was  successfully  applied  to  a  system  of 
four  nonlinear  partial  differential  equations  (mass  conservation,  two  components  of  linear 
momentum  equation,  and  energy  equation)  governing  unsteady  two-dimensional  flow  of 
compressible,  rotational,  inviscid  fluid.  This  system  is  known  as  Euler  equations  of  gas 
dynamics.  The  basic  integration  algorithm  was  an  explicit  scheme  that  utilizes  Runge- 
Kutta  time-stepping  and  finite  volume  formulation  for  spatial  discretization.  The 
algorithm  is  known  as  Jameson’s  scheme  and  represents  the  fastest  presently  available 
integration  method  for  Euler  equations  of  gas  dynamics. 

A  new  Distributed  Minimal  Residual  (DMR)  method  for  the  acceleration  of  explicit 
iterative  algorithms  for  the  numerical  solution  of  systems  of  partial  differential  equations 
has  been  developed  by  Lee  and  Dulikravich.  The  method  is  based  on  the  idea  of  allowing 
each  partial  differential  equation  in  the  system  to  approach  the  converged  solution  at  its 
own  optimal  speed  while  at  the  same  time  communicating  with  the  rest  of  the  equations 
in  the  system.  The  DMR  method  belongs  to  a  general  class  of  the  extrapolation 
techniques  in  which  the  solution  is  updated  using  information  from  a  number  of 
consecutive  time  steps  in  such  a  way  that  the  L2  norm  of  future  residual  is  minimized. 
Unlike  in  other  similar  methods,  each  component  of  the  solution  vector  is  updated  using 
a  separate  sequence  of  acceleration  factors.  The  idea  of  using  different  acceleration 
factors  for  each  component  of  a  solution  vector  is  similar  to  that  of  dynamic 
preconditioning.  This  allows  each  equation  to  evolve  at  its  own  optimal  convergence 
rate.  Moreover,  the  acceleration  factors  are  determined  from  the  governing  equations  so 
that  only  a  few  consecutive  solutions  are  required  for  a  successful  application  of  the 
DMR  method.  This  acceleration  scheme  was  applied  to  the  system  of  time-dependent 
Euler  equations  of  inviscid  gasdynamics  in  conjunction  with  the  finite  volume  Runge- 
Kutta  explicit  time-stepping  algorithm.  Using  DMR  without  multigridding,  between  30% 
and  70%  of  the  total  computational  efforts  were  saved  in  the  subsonic  compressible  flow 
calculations.  The  DMR  method  seems  to  be  especially  suitable  for  stiff  systems  of 
equations  and  can  be  applied  to  other  systems  of  differential  equations  and  other 
numerical  algorithms.  Specifically,  the  DMR  method  was  applied  to  an  artificial 
compressibility,  explicit,  Runge-Kutta  time  stepping  algorithm  for  steady,  incompressible, 
Navicr-Stokes  equations.  A  two-dimensional  analysis  computer  code  in  a  generalized 
curvilinear  coordinate  system  was  developed  and  its  accuracy  has  been  compared  to 
known  numerical  solutions.  The  algorithm  has  been  successfully  accelerated  using  the 
DMR  method,  resulting  in  25%-70%  reduction  in  computing  time. 
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II.  STATUS  OF  THE  RESEARCH 


The  objective  of  the  research  project  was  to  provide  a  sound  mathematical  theory 
for  non-linear  iterative  acceleration  schemes  using  multiple  optimal  acceleration  factors 
and  to  test  the  method  on  several  non-linear  differential  systems.  Particular  emphasis 
was  placed  on  developing  computer  programs  for  accelerating  the  convergence  and 
enhancing  stability  of  iterative  solutions  of  the  non-linear  systems  of  partial  differential 
equations  of  fluid  mechanics. 

During  the  course  of  this  research  project,  both  analytical  and  software  development 
aspects  were  addressed.  A  general  theory  of  optimal  acceleration  factors  for  the  multi- 
step  iterative  solution  of  systems  of  non-linear  partial  differential  equations  based  on  the 
minimal  residual  concept  was  developed  with  the  special  emphasis  on  mixed-type  fvst  ms 
of  partial  differential  equations.  The  new  methods  was  tested  on  a  variety  of  practical 
examples  governed  by  the  compressible  two-dimensional  inviscid  flow  equations  (Euler 
equations)  and  viscous  incompressible  laminar  flow  equations  (Navier-Stokes  equations). 
Subsonic  and  transonic  flow  fields  were  calculated  for  geometries  including  nozzles, 
airfoil  cascades,  airfoil  in  an  unbounded  domain,  and  a  driven  cavity  problem. 

Two  graduate  students,  Mr.  Chung-Yuan  Huang  and  Mr.  Stephen  R.  Kennon,  have 
finished  their  doctorate  degrees  in  the  summer  of  1987  at  the  University  of  Texas  at 
Austin,  while  supervised  by  Professor  George  S.  Dulikravich  from  Penn  State  University 
who  continued  to  serve  as  an  adjunct  faculty  with  the  University  of  Texas.  Both  Dr. 
Huang  and  Dr.  Kennon  were  partially  supported  by  this  grant  and  the  preceding  grant 
from  the  AFOSR/NM  with  Professor  David  M.  Young  as  co-principal  investigator.  This 
fact  was  acknowledged  in  their  doctoral  dissertations. 

Dr.  Huang  has  applied  GNLMR  to  a  number  of  scalar  nonlinear  partial  differential 
equations  and  to  a  system  of  Euler  equations  of  gasdynamics.  He  now  works  as  a 
postdoctoral  Research  Scientist  with  Professor  J.  Tinsley  Oden  at  the  University  of 
Texas.  Dr.  Kennon  has  developed  a  number  of  new  ideas  for  finite  elements  in 
gasdynamics  including  acceleration  of  iterative  algorithms.  He  now  works  as  an  Assistant 
Professor  in  the  Aerospace  Engineering  Department  of  the  University  of  Texas  at 
Arlington. 

A  new  graduate  research  assistant  at  Penn  State  was  supported  with  the  grant 
AFOSR-87-OI21.  Mr.  Seungsoo  Lee  is  a  Ph.D.  candidate  who  has  developed  the  DMR 
method  and  applied  it  to  the  explicit  finite  volume  Runge-Kutta  scheme  (Jameson’s 
algorithm)  for  Euler  equations  of  gasdynamics.  He  has  derived  all  the  governing 
equations  in  a  fully  conservative  nondimensionalized  form  suitable  for  discretization  on  a 
general  nonorthogonal  curvilinear  computational  grid. 

Recently,  Mr.  Lee  has  successfully  implemented  the  DMR  concept  in  an  explicit 
algorithm  for  the  numerical  integration  of  Navier-Stokes  equations  of  laminar, 
incompressible  flows  through  nozzles  and  cascades. 

Another  Ph.D.  candidate,  Mr.  Daniel  J.  Dorney,  worked  on  the  analysis  of  existing 
numerical  dissipation  models  and  on  physically  based  dissipation  formulations  for  Euler 
and  Navier-Stokes  equations  of  gasdynamics.  Mr.  Lee  and  Mr.  Dorney  worked  together 
on  implementing  DMR  in  both  Euler  and  Navier-Stokes  codes  with  the  physically  based 
artificial  dissipation. 

A  Visiting  Research  Scientist,  Mr.  Ren  Bing,  was  involved  on  the  project  at  the 
Penn  State  Uni  ersity  for  six  months.  He  performed  a  thorough  survey  of  all  Total 
Variation  Diminishing  (TVD)  type  schemes  for  controlling  numerical  dissipation. 
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Dynamics."  Speaker:  Prof.  G.  S.  Dulikravich. 

2.  Dulikravich,  G.  S.,  "Optimization  of  Explicit  Multi-Step  Algorithms,"  paper 
presented  at  the  First  International  Conference  on  Industrial  and  Applied 
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of  different  artificial  dissipation  models.  Prof.  G.  S.  Dulikravich. 
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G.  S.  Dulikravich. 
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discoveries,  inventions  or  patents  resulting  from  this  research  project. 

VIII.  SUGGESTIONS  FOR  FUTURE  RESEARCH 

We  also  developed  a  DMR  version  of  a  transonic  Navier-Stokes  finite  volume  code 
for  two-dimensional  shock/boundary  layer  interaction.  In  addition,  Mr.  Lee  implemented 
DMR  in  a  two-dimensional  implicit  ADI  (Beam-Warming)  solver  for  incompressible  Navier- 
Stokes  equations.  He  also  developed  a  fully  three  dimensional  DMR  version  of  an  explicit 
code  for  incompressible  Navier-Stokes  equations.  These  three  codes  remain  to  be  tested 
especially  when  using  DMR  formulation  with  implicit  algorithms. 

Notice  that  all  numerical  results  with  the  DMR  method  were  obtained  without  the 
standard  acceleration  techniques  such  as  explicit  and  implicit  residual  smoothing,  enthalpy 
damping,  multigridding  and  vectorization.  These  methods  could  be  combined  with  the 
DMR  to  further  accelerate  the  iterative  algorithms. 

Furthermore,  it  would  be  highly  desirable  to  study  the  effect  of  grid  clustering,  grid 
size,  grid  orthogonality  and  grid  structure  on  the  DMR.  In  addition,  domain  partitioning 
could  be  used  with  different  DMR  sequences  in  each  domain  thus  leading  to  accelerated 
parallel  processing  capabilities. 
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The  Generalized  Nonlinear  Minimal  Residual  (GNLMR)  method  is  shown  to  consistently  acceler¬ 
ate  and  stabilize  iterative  algorithms  for  solving  nonlinear  problems  by  using  the  optimized  explicit 
multistepping.  The  examples  presented  in  this  paper  illustrate  the  beneficial  effects  of  the  optimized 
multistep  algorithm  on  the  computational  efficiency  and  the  convergence  rate  as  applied  to  several 
nonlinear  problems  in  fluid  dynamics.  The  significant  reduction  in  computing  lime  when  using  the 
multiple  optimized  acceleration  factors  is  only  negligibly  weighed  down  by  the  computation  costs  due 
to  the  requirements  for  additional  computer  storage. 


1.  Introduction 

The  relaxation  factor  used  in  accelerating  an  iterative  method  to  obtain  the  converged 
solution  plays  the  same  role  as  the  time  step  size  in  advancing  the  transient  solution  to  the 
steady-state  solution  for  a  time-dependent  problem.  The  classical  analyses  for  the  stability  of 
numerical  schemes  for  solving  time-dependent  problems  neglect  boundary  conditions  and 
assume  a  uniform  computational  grid.  Furthermore,  these  analyses  are  based  on  linear 
equations  with  constant  coefficients  and  the  assumptions  of  small  perturbations  and  ap¬ 
plicability  of  Fourier  analysis  [1,2].  However,  Cheng  [2]  pointed  out  that  the  perturbation  of 
the  error  in  the  finite  difference  calculations  may  not  be  small  and  that  the  error  in  the  finite 
difference  calculations  may  not  satisfy  the  conditions  for  Fourier  series  expansion.  In  addition. 
Mitchell  and  Griffiths  [1]  pointed  out  that  the  errors  due  to  approximate  or  additional 
boundary  conditions  are  represented  by  modes  which  are  not  of  Fourier  type.  Thus,  the  linear 
stability  analysis  usually  results  in  overly  restrictive  and  even  incorrect  conclusions. 

The  numerical  experiments  performed  by  Kennon  and  Dulikravich  [3]  and  Kennon  [4] 
using  the  NonLinear  Minimal  Residual  (NLMR)  method  showed  that  the  usual  Courant- 
Friedrichs-Lewy  (CFL)  number  limitation  for  both  linear  and  nonlinear  problems  can  be 
significantly  exceeded.  The  NLMR  method  provided  a  simple  analytic  way  to  determine  the 
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optimal  acceleration  factors  for  both  linear  and  nonlinear  problems.  However,  the  elementary 
time  steps  used  for  obtaining  the  corrections  still  follow  the  CFL  number  limitation  concluded 
from  the  linear  analysis. 

The  generalized  nonlinear  minimal  residual  (GNLMR)  method  developed  by  Huang, 
Kennon  and  Dulikravich  [5]  provided  a  practical  analytical  tool  to  determine  the  exact 
stability  conditions  for  both  linear  and  nonlinear  problems  in  arbitrary  domains.  If  accurate 
time  evolution  is  required  when  solving  an  unsteady  problem,  the  limitation  on  the  time  step 
size  can  be  analytically  determined  by  using  the  GNLMR  method.  If  transient  behavior  is  of 
no  interest,  the  GNLMR  method  can  be  applied  to  determine  the  optimal  value  of  the  time 
step  size  (optimal  acceleration  factor)  to  minimize  the  number  of  time  steps  (number  of 
iterations)  for  obtaining  the  steady-state  converged  solution. 

The  main  objective  of  this  paper  is  to  investigate  the  effects  of  the  optimized  n  itistep 
algorithm  on  the  computational  efficiency  and  on  the  monotonicity  of  convergence  rate  of  the 
GNLMR  method.  The  analytic  investigation  is  confirmed  on  four  nonlinear  test  cases:  the 
one-dimensional  and  two-dimensional  viscous  Burgers’  equations  and  the  two-dimensional 
incompressible  and  compressible  Stream-Function-Coordinate  (SFC)  equations  [6], 

2.  Theoretical  aspects 

2.1.  Multistep  minimal  residual  method  for  linear  problems 
Let  us  first  consider  a  well-posed  linear  initial  value  problem: 

dtp! dr  =  L<p  -  F  in  D  , 
tp  =  ipB  on  dPi  , 

<P  =  <Po  at  r  =  r0 

Define 

r'  =  l<p'  —  f 

as  the  residual  vector  at  time  level  r.  Here,  l  denotes  the  scheme-dependent  discrete  analog  of 
L,  /  is  the  discrete  analog  of  F  and  also  includes  boundary  terms. 

Assume  that  M  steps  are  used  to  iterate  at  each  time  level  t.  Using  the  Einstein  summation 
convention  where  repeated  subscripts  are  summed,  the  multistep  algorithm  for  (1)  is  then 
defined  as  follows: 

<p,+‘  =  <p‘  +  a)m8m  ,  m  =  l,2,  (3) 

where 

5i  *  V-/, 

5m  =  r-‘(51),  m>  1,  (4) 

are  the  corrections  at  step  m.  Coefficients  a>m  are  the  corresponding  relaxation  factors  to  be 
determined  by  minimizing  the  L2  norm  of  the  residual  at  time  level  (t  +  1).  With  the  definition 


(1) 

(2) 
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of  residual  vector  (2)  and  the  help  of  (3),  the  following  relation  can  easily  be  verified: 

f-t  1  t  *  fc  /  .  tm  i 

r  =  r  +  coj 8m  =  r  +  u)J  r  .  V0) 

The  L;  norm  of  the  residual  vector  at  time  level  (t  +  1)  is  then 

Ik" 'll'  =  jlr'H'  +  2wm(r',  I8m)  +  ( 18  m .  I8n)wmwn  ,  m,  n  =  1.2 . A/  .  (6) 

It  should  be  pointed  out  that  the  boundary  condition  for  the  corrections  8m  in  (3)  is 
apparently  zero.  However,  the  boundary  conditions  ft  .  the  residua!  vector  r  and  corrections 
8m  in  (5)  can  be  determined  either  b.  extrapolation  from  the  interior  points  or  simply  bv 
setting  them  equal  to  zero.  The  resiciual  norm  will  then  converge  to  the  norm  of  the  truncation 
error  of  the  difference  scheme  if  the  first  method  is  applied  and  to  the  machine  accuracy  if  the 
second  method  is  applied. 

The  highest  rate  of  convergence  is  possible  when  wm  are  the  solutions  of  the  following 


system  of  linear  equations: 

ar/awm=  0  or  (r',l8m)  +  (l8m,l8n)wn  =  Q,  (7) 

where  the  rate  of  convergence  P  is  defined  as 

r=-log(||r'  +  1||/||r'||).  (8) 

Multiplying  (7)  by  wm,  it  follows  that 

(r‘,l5„)(om  +  (l8m,I8n)umwn  =  0 .  (9) 

Subtracting  (9)  from  (6)  and  using  (7)  results  in 

Ill'll2  -  Ik'll2  =  (r\  I8m)wm  =  - (I8m .  l8n)o>m»n  =  -  ( dQ  <0  .  (10) 

Jit 


Thus,  the  residual  norms  for  the  multistep  minimum  residual  method  show  a  monotone 
convergence  behavior  which  guarantees  the  stability  of  the  iterative  scheme  and  produces  the 
highest  rate  of  its  convergence. 

2.2.  Optimization  of  the  Euler  scheme  for  nonlinear  problems 

For  clarity,  we  consider  two-dimensional  problems  and  c  quations  in  conservative  form  only. 
The  extension  to  multidimensional  problems  and  nonconservative  equations  is  then 
straightforward. 

rin-  conservative  fonn  of  the  governing  equations  for  most  engineering  problems  can  be 
written  as: 

dcpldr  =  LlfN"((p,  tpx,  <pv)  -  F  ,  ('ll) 
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where  the  operators  are 

L,  =  bl bx  ,  L,  =  d/by  , 

and  A'1'  is  the  nonlinear  differential  operator  in  coordinates  .v, .  Using  the  Euler  one-step, 
time-consistent,  explicit  scheme,  the  finite  difference  form  of  (tl)  can  be  written  as: 


<f  "  1  =  ip'  +  A rr‘ ,  (12) 

where 

r  =  if[,  <p\ )  —  f  (13) 

is  defined  as  the  residual  at  time  level  t.  Therefore,  the  residual  at  time  level  (t  +  1)  can  be 
expressed  as 

r  =  /,.  A  (<p  .  tfx  ,  <pv  )  -  f  .  (14) 

After  expanding  the  nonlinear  discretized  operator  N"  in  a  Taylor  scries,  it  follows  that 

r'M  =  UA'V.  if‘x.  <p\)  +  [(bN‘/d<p'))r:  +  (bN‘ib<p'x){r')x 

+  (bN‘7 bip[  )(r')vjAr  +  0(A72)}  -  /.  (15) 

In  summary, 

r'*1  =  r'  +  a/,(Ar)p,  1  ^  p  S  P  ,  (16) 


where  P  is  the  degree  of  the  nonlinearity  of  the  operator  N.  Equation  (16)  indicates  that  the 
res:dual  at  time  level  (f  +  1)  is  a  polynomial  (henceforth  called  Residual  Polynomial  [3]  or  RP) 
of  the  time  step  size,  Ar.  Thus,  the  L2  norm  of  the  residual  at  time  level  (t  +  1)  can  be 
expressed  as 

||r'Tl||:  =  ||r'||:  +  2(r'ap)(Arr +  (n/,.fl,)(A7r(ATr  .  l*p.  q  *£  P  ■  (17) 

Equation  (17)  implies  that  the  residual  norm  at  time  level  (r  +  1)  is  a  positive  polynomial 
(henceforth  called  Minimizing  Polynomial  [3]  or  MP)  of  the  time  step  size  Ar,  which  is  to  be 
determined.  Thus,  the  convergence  of  scheme  (12)  will  be  guaranteed  provided  that  Ar  is 
chosen  in  such  a  way  that  T  >  0.  The  highest  rate  of  convergence  can  be  achieved  only  when 
At  is  chosen  as  the  optimizer  of  the  minimizing  polynomial  (17)  such  that  |j  rr * 1 1(  is  an 
infimum.  However,  the  determination  of  the  optimizer  needs  special  numerical  techniques  [7], 
To  avoid  this  difficulty,  the  linearized  operator  (3-5]  of  Nr  may  be  applied.  If  Nl  is  truncated 
to  the  first  order  in  Ar  (linearized  operator),  the  approximate  residual  vector  is 

r*'*1  -  r‘  +  a,Ar  ,  (18) 

a,  =  U(dtf7a? V  +  O.V7a^)(r')x  +  (0A'7d<)(r'),]  .  (19) 


where 
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Then,  the  approximate  MP  is 

Ik*'* 'll’  =  Ik'll'  +  2(a,.  r')Ar  +  (a,,  fl.KAr)2  .  (20) 

The  optimal  time  step  for  the  explicit  Euler  scheme  can  be  easily  found  as 

(At)opi  =  r')/||a,||;  .  (21) 

2.3.  The  generalized  nonlinear  minimum  residual  ( GNLMR )  method 

The  GNLMR  method  actually  is  the  application  of  the  methods  described  in  the  previous 
sections.  The  multistep  algorithm  for  nonlinear  problems  is  defined  as 

=  <p‘  +  wm8m  +  0(w^,)  .  m  =  1.2 . M.  (22) 

where  repeated  indices  are  summed.  The  correction  at  the  first  step  if  defined  as 

8,  =  r'  =  /..ATVVtVv)-/.  (23) 

The  correction  at  step  m  >  1  is  defined  as 

8.  =  l,.[(dNv/d'P,)8m_i  +(3AT/V,  )(«*-,),  +  (dAT/d*>,J(8m_1)v]  •  (24) 


The  coefficients  of  the  higher-order  terms  of  a)m  can  be  obtained  by  Tavlor-series  expansion.  If 
only  linear  terms  of  wm  are  retained,  the  residual  polynomial  (RP)  at  time  level  (t  +  1)  can  be 
expressed  by  Taylor-series  expansion  as 

=  +  <p\  +  wm(8m)t,  tp\.  +  wm(5m)J  -/ 

=  r'  +  ^{[(8JV7ay')6m  +  (dN"ldtp[)(8m)x  +  (dNvlhip\ )(8^)Jwm  +  Ofw;)}  . 

(25) 

Therefore,  the  minimizing  polynomial  (MP)  at  time  level  (t  +  1)  can  be  determined  as 

IK'll:  =  lkT  +  s(<o.  (26) 

where  g(wm)  is  a  polynomial  in  (um.  For  a  highly  nonlinear  differential  equation,  g  will  be  a 
complicated  multivariable  polynomial  that  depends  on  the  total  number  of  intermediate  steps 
M  that  were  used  and  the  degree  of  the  nonlinearity  of  the  differential  operator  N".  Thus,  a 
fast  and  accurate  procedure  of  determining  the  optimizer  of  MP  is  required  for  the  GNLMR 
method  to  guarantee  the  highest  rate  of  convergence.  If  the  linearized  operator  of  N"  is  used, 
the  method  that  was  described  in  Section  2. 1  can  be  applied  to  determine  the  approximate 
optimizer  of  (26). 

The  GNLMR  method  requires  (M  +  1)  times  larger  computer  storage  to  save  the  correc- 
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tions  from  M  intermediate  steps  than  does  the  single-step  nonaccelerated  scheme.  Some 
additional  algebraic  operations  are  also  required  to  determine  the  coefficients  of  the  MP  which 
are  obtained  by  integrating  the  corrections  over  the  entire  domain.  Th°  storage  requirement 
of  the  GNLMR  method  is  quite  acceptable  when  compared  with  the  excessive  storage 
required  by  the  GMRES  method  [8.9],  It  should  be  pointed  out  that  the  GMRES  method 
also  needs  a  large  number  of  arithmetic  operations  not  only  for  orthonormalizing  search 
directions  but  also  for  determining  the  optimal  weighing  parameters  in  updating  the  iterative 
solutions. 


3.  Numerical  examples 

Four  test  cases  were  used  to  demonstrate  the  application,  the  computational  efficiency,  and 
the  monotone  convergence  behavior  of  the  GNLMR  method.  Since  it  was  found  [5]  that  the 
linearized  residual  polynomial  still  guarantees  a  relatively  high  convergence  rate,  it  will  be 
used  for  all  test  cases.  The  first  two  cases  representing  the  one-dimensional  and  two- 
dimensional  viscous  Burgers'  equation  were  solved  by  the  time-dependent  technique  as 
described  in  Section  2. 

The  last  two  test  cases,  the  two-dimensional  incompressible  and  compressible  stream- 
function-coordinate  (SFC)  equations  [6]  were  solved  in  their  steady-state  and  nonconservative 
form.  Liebman's  or  Gauss-Seidel's  method  was  applied  to  determine  the  correction  at  each 
intermediate  iterative  step  m. 

Details  about  the  control  parameters  such  as  grid  size,  stopping  criteria,  and  number  of 
acceleration  factors  used  for  each  test  case  are  summarized  in  Table  1.  For  all  test  cases, 
comparisons  are  based  on  the  relative  improvement  of  computational  efficiency  that  can  be 


Table  1 


Summary  of  the  control  parameters  for  numencal  test  cases 


Test  case 

Max. 
no  of 

(a  used 

Boundary  conditions 
for  residual  and 
corrections 

Stopping 

criteria 

Grid 

size 

Case  1: 

ID  Burgers' 
equation 

8 

Extrapolation  from 
interior  data 

Zero 

r  « 10  " 

l|r'||«10  * 

41 

Case  2: 

2D  Burgers’ 
equation 

S 

Extrapolation  from 
interior  data 

Zero 

r«  io  * 

llr-ll*  10-* 

51  X  51 

Case  3: 

2D  incomp. 

SFC  equation 

8 

Zero 

Hr'H  «  III  " 

47  x  II 

Case  4: 

2D  comp. 

SFC  equation 

8 

Zero 

||r'H  «  10  " 

61  x  11 

10 
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obtained  using  different  total  number  of  intermediate  steps  M.  The  relative  improvement  of 
computational  efficiency  is  defined  as 

r,u  =  T0iTM.  (27) 

where  Tn  denotes  the  computing  time  spent  for  the  nonaccelerated  method  and  7\,  denotes 
the  computing  time  spent  for  the  acceleration  method  based  on  the  stopping  criteria  as 
described  in  Table  1.  The  results  are  summarized  in  Fig.  9. 

5. 1 .  Burgers '  equations 

According  to  the  notations  defined  in  Section  2,  the  one-dimensional,  viscous  Burgers' 
equation  can  be  written  as 

=  a/ajc[7V(sP,  «pr)j  -  (28) 

where 

N(<p.  <f()  =  -  ;<p'  +  vip,  .  (29) 

and  v  is  the  viscosity  coefficient.  In  this  example  v  =  0.07  is  used.  The  initial  and  the  boundary 
conditions  are  chosen  as  follows: 

<p(l,r)  =  0.  <p(0.  r)  =  1  .  q>{x.  0)  =  1  -  .t  .  (30) 

The  two-dimensional,  nonlinear,  viscous  Burgers'  equation  solved  by  Ghia  et  al.  [10]  was 
chosen  in  the  presented  test  case.  According  to  the  notations  defined  in  Section  2,  it  can  be 
expressed  as 

du/dr  =  d/dx[Nl(u.  uj]  +  d/dy[N:(u.  wv)]  .  (31) 

where 

/Vl(w,  uz)  =  ux  -  Ay(iw;  -  uV) ,  N2(u.  u})  =  u\  -  \x(\u2  -  uil) .  (32) 

Here,  A  is  a  parameter  and  U  is  a  constant.  The  values  of  U  and  A  used  in  this  test  case  are  0.5 
and  2.0,  respectively. 

The  FTCS  scheme  is  applied  to  discretize  (28)  and  (31).  If  the  linearized  form  of  operator 
N  is  used  with  M  steps  at  each  time  level  t.  the  residual  polynomial  is  truncated  up  to  its  first 
order  as 

RP  =  r'*1  =  r'  +  amwm  ,  (33) 

where 

flm  =  a/d*[-«p'5m  +  Kf>J,]  (34) 

for  the  one-dimensional  case  and 

=  d/dx[(dNl/du)8m  +  (a/V'/dut)(8m)J  +  dldy[(dN2ldu)8m  +  (dN2/  a«v)(8m)v] 

(35) 


22 


G.S.  Dulikravich,  C.-Y.  Huang,  Fast  algorithms  based  on  time  stepping 


for  the  two-dimensional  case.  The  corrections  5„,  at  each  intermediate  step  m  can  be 
determined  by  (23),  (24).  The  minimizing  polynomial  MP  for  both  cases  is  then 


MP®  ||r'*1||J  =  ||r'||-  +  2(r'.  am)wm  +  (am,  an)wmwn  . 


Thus,  the  optimal  acceleration  parameters  tom  can  be  easily  determined  by  solving  the 
following  system  of  linear  equations: 


A  at  =  b  .  (37) 

where 

=  (om.an) .  (38) 


bm  =  ~(r\  aj 


and  Amn  is  a  symmetric  matrix  of  order  M. 

The  boundary  conditions  of  the  residual  and  corrections  in  (33)  were  determined  by 
extrapolating  them  from  the  interior  points  (Case  a)  or  by  explicitly  enforcing  them  to  be  zero 
(Case  b).  The  stopping  criteria  for  all  cases  were  such  that  the  computations  were  terminated 
when  the  asymptotic  rate  of  convergence  (Case  a)  or  the  norm  of  the  residual  (Case  b) 
approached  the  machine  accuracy. 

It  must  be  mentioned  that  the  norm  of  the  truncation  error  represents  the  maximum 
attainable  accuracy  of  a  numerical  scheme  and  is  obviously  scheme-dependent.  It  can  be  seen 
from  Figs.  2(a)  and  4(a)  that  under  the  same  stopping  criteria  the  residual  norms  for  all  cases 
converge  to  the  values  corresponding  to  the  respective  norms  of  the  truncation  error. 
Moreover,  Figs.  1(a)  and  2(a)  illustrate  that  the  accuracy  of  the  nonaccelerated  scheme  can  be 
improved  by  applying  the  GNLMR  method. 

Since  the  linearized  operators  were  used  in  these  two  test  cases,  the  convergence  history 
shown  in  Figs.  1-4  exhibit  a  similar  behavior  as  in  the  linear  problems  as  solved  in  our  earlier 
works  [5].  It  is  obvious  that  if  the  GNLMR  method  is  applied,  the  number  of  iterations  and 
the  computing  time  required  to  achieve  the  asymptotic  rate  of  convergence  are  considerably 
lowered  as  compared  to  the  nonaccelerated  schemes  (Figs.  1(a),  2(a),  3(a)  and  4(a)). 
Moreover,  the  time  required  by  the  GNLMR  method  for  marching  the  solution  from  the 
asymptotic  state  to  a  fully  converged  solution  is  much  shorter  than  with  the  nonaccelerated 
method.  The  improvement  of  the  computational  efficiency  that  can  be  obtained  using  a 
different  number  of  intermediate  steps  M  is  summarized  in  Fig.  9. 

Although  the  computational  efficiency  increases  significantly  with  the  increasing  number  of 
intermediate  steps  A/,  the  improvement  becomes  less  pronounced  and  even  shows  a  reverse 
trend  after  approximately  M  =  5  in  the  one-dimensional  case.  The  reason  for  this  unexpected 
result  is  that  when  using  a  multistep  algorithm,  an  M  x  M  matrix  has  to  be  inverted  (directly) 
at  each  time  level,  t.  The  number  of  operations  and  computing  time  required  for  the  direct 
inversion  of  a  matrix  grows  very  fast  with  the  increase  of  the  matrix  size,  thus  countering  the 
benefits  of  adding  more  intermediate  steps  in  the  multistep  procedure  especially  for  one¬ 
dimensional  problems. 
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L OG 1 0 ( R£ SNOR  M)  VS.  N(ITER)  ,  1-0  BURGERS'  EQUATION 
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Fig.  1(a).  Residual  norm  versus  the  number  of  iterations.  ID  Burgers'  equation. 


3.2.  Stream-function-coordinate  (SFC)  equations 

The  two-dimensional  stream-function-coordinate  (SFC)  equation  for  an  irrotational,  in  vis¬ 
cid,  steady  flow  derived  by  Huang  and  Dulikravich  [6]  is  given  by: 

( yl  ~  eh,*  ~  2y,y*y,#  +  (1  +  y])yM  =  0  ,  (40) 

where  cr  represents  the  compressibility  and  is  equal  to  zero  for  incompressible  flows.  It  is 
defined  as 

o  =  {p*a*  Ipaf  =  \\(y  +  l)-*(y-l , 


(41) 
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Fig.  1(b).  Residual  norm  versus  ihe  number  of  iterations.  ID  Burgers'  equation. 


where  p  and  a  denote  the  local  density  and  the  local  speed  of  sound,  respectively,  and  y  is  the 
ratio  of  the  specific  heats.  The  superscript  terms  denote  the  characteristic  quantities  of  the 
flow.  It  can  be  shown  [6]  that  cr  is  an  implicit  function  of  y(  and  y ii,  that  is 

(i  +y;)'yi  =  [(y  +  (y+n  -  2]/[(y  -  1M .  (42) 

Let  us  define 

ci  =  y,  ’  c:  =  y *,  c3  =  yxx, 

'  cs  =  v„*  •  (43) 


// 


Then  (40)  can  be  rewritten  as 
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Fig.  2(a).  Residual  norm  versus  the  computing  time.  ID  Burgers'  equation 

iV(c, ,  c,,  c3,  c4.  c5 )  =  0  ,  (44) 

or 

[C;  -  cr(c,,  c,)]c,  -  2ctc: c\  +  ( l  +  c()c<  =  0  .  (45) 

Assume  that  a  uniform  computational  grid  is  used  (both  A.v  and  At/r  are  constant)  and 
central  differencing  is  applied  to  discretize  all  derivatives.  The  finite  difference  approximation 
of  the  SFC  equation  can  be  expressed  as 

y,.,  =  [(>'*  ~  er)/AjT  +  ( 1  +  >;)/At/r:]"'[(y;  -  cr)( +  V, _ ,  .,)/(2A.r:) 

+  d  +  .Vr)(  y,.,.,  +  y, .,-,)/(2At/r;)  -y,ytyxt]  ■ 


(46) 


where  t  represents  the  iteration  level,  w  is  the  relaxation  factor,  and  ;  is  the  correction  at 
iteration  level  /.  It  is  defined  as 
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Fig.  3(a).  Residual  norm  versus  the  number  of  iterations.  3D  Burgers'  equation. 


Here,  y {*'  is  the  temporary  value  of  y  at  iteration  level  (t  +  1)  obtained  by  applying  (46)  with 
any  fundamental  iterative  scheme.  In  the  presented  studies,  Liebman's  method  (a>  =  1)  was 
used  as  the  fundamental  iterative  scheme  and  will  be  henceforth  referred  to  as  the  nonacceler¬ 
ated  method.  Assume  that  M  steps  are  used  in  the  GNLMR  method.  The  solution  is  then 
updated  by  using 


<+/  =</  +  ,  m  =  l,2 . M  , 


(49) 


n 
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Fig.  3(b).  Residual  norm  versus  the  number  of  iterations,  2D  Burgers'  equation. 


where  8m  are  the  corrections  at  each  intermediate  step  m.  They  are  obtained  by  successively 
applying  Liebman's  method.  The  optimal  values  of  u>m  based  on  a  linearized  RP  can  be 
determined  by  solving  (37)  with 


r'  =  N(ct,  c2,  Cj,  ct,  c5)r , 


(50) 


=  [(dN/dCj)(8m)x  +  (dN/dc2)(8m),  +  ( dN/dc3)(8Jxil , 
+  (alV/dc4)(6m)„  +  (dN/dcs)(SJ„]. 


(51) 


4 


For  the  incompressible  case,  a  uniform  flow  around  a  cascade  of  doublets  was  solved  [6], 
For  the  compressible  case,  a  subsonic  flow  with  free-stream  Mach  number  ,M,  =  0.65  around  a 
NACA  0012  airfoil  in  a  channel  with  height/chord  ratio  =3.6  was  solved  [6], 

Since  linearized  operators  were  used  in  these  two  cases,  and  the  boundary  conditions  for 
the  residual  and  corrections  in  (33)  were  set  to  zero,  the  residual  norm  will  converge  to 
machine  accuracy.  Therefore,  the  stopping  criteria  for  these  two  cases  was  chosen  in  such  a 
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Fig.  4(b).  Residual  norm  versus  the  computing  time,  2D  Burgers'  equation. 


way  that  as  the  residual  norm  approaches  machine  accuracy,  the  computation  is  forced  to 
stop. 

The  improvement  of  the  computational  efficiency  is  summarized  in  Fig.  9.  Both  cases  show 
that  the  computational  efficiency  is  increased  significantly  by  increasing  the  number  of 
intermediate  steps  M. 

The  numerical  results  for  these  two  cases  are  summarized  in  Figs.  5-8.  Both  cases  exhibit  a 
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common  feature  that  the  monotonicitv  and  the  rate  of  convergence  increase  with  the  increase 
of  the  total  number  of  intermediate  steps  .V/.  Most  importantly.  Figs.  6  and  8  show  that  with  a 
specified  minimum  computing  time,  the  difference  in  the  residual  norms  between  the 
nonaccelerated  method  and  the  GNLMR  method  varies  between  one  to  eight  orders  of 
magnitude  depending  on  the  number  of  steps  used  in  the  GNLMR  method.  This  fact  strongly 
proves  the  computational  efficiency  that  can  be  obtained  using  the  GNLMR  method  [8). 


4.  Concluding  remarks 

Four  numerical  test  cases  for  nonlinear  problems  in  fluid  dynamics  were  presented  to 
demonstrate  the  applicability,  computational  efficiency,  and  monotone  convergence  behavior 
of  the  GNLMR  method.  It  was  found  that  even  though  the  theory  of  the  GNLMR  method  is 
based  on  the  evolution  problems  and  equations  in  conservative  form,  the  method  can  be 
applied  equally  successfully  to  the  solutions  of  steady-state  problems  governed  by  equations  in 
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Fig.  7.  Residual  norm  versus  the  number  of  iterations,  compressible  SFC  equation. 


nonconservative  form.  The  results  for  all  test  cases  show  that  when  applying  the  GNLMR 
method  to  nonlinear  problems,  the  number  of  iterations  and  the  corresponding  computer  time 
are  considerably  lowered  by  increasing  the  number  of  intermediate  time  steps. 

Since  the  explicit  multistep  algorithm  was  employed  in  developing  the  GNLMR  method, 
the  advantage  of  accelerating  the  convergence  rate  of  the  iterative  process  is  partially  offset  by 
some  extra  costs.  These  are  caused  by  the  requirements  for  additional  storage  in  order  to  save 
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Fig.  8.  Residual  :;orm  versus  the  computing  time,  compressible  SFC  equation. 


corrections  obtained  from  each  intermediate  step  and  by  the  additional  arithmetic  operations 
to  determine  the  coefficients  of  the  minimizing  polynomial.  In  practice,  a  maximum  gain  in 
computational  efficiency  can  be  obtained  with  a  moderate  number  (usually  not  more  than  five) 
of  intermediate  steps.  The  requirement  for  additional  storage  linearly  increases  with  the 
number  of  intermediate  time  steps  used  and  represents  only  a  fraction  of  the  computer  storage 
required  by  the  GMRES  method  [8]. 
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A  new  algorithm  for  the  acceleration  of  explicit  iterative  schemes  for 
the  numerical  solution  of  nonlinear  systems  of  partial  differential 
equations  has  been  developed.  The  method  is  based  on  the  idea  of  allowing 
each  partial  differential  equation  in  the  system  to  approach  the  converged 
solution  at  its  own  optimal  speed.  The  DMR  (Distributed  Minimal  Residual) 
method  introduces  a  separate  sequence  of  optimal  weighting  factors  to  be 
used  for  each  component  of  the  general  solution  vector.  The  acceleration 
scheme  was  applied  to  a  highly  nonlinear  coupled  system  of  four 
time-dependent  partial  differential  equations  of  inviscid  gasdynamics  in 
conjunction  with  the  finite  volume  Runge-Kutta  explicit  time-stepping 
algorithm.  Using  DMR  without  mu  1 1 igr idding ,  between  30%  and  70%  of  the 
total  computational  efforts  were  saved  in  the  subsonic  compressible  flow 
calculations.  DMR  method  offers  most  time  savings  when  applied  to  stiff 
systems  of  equations. 

Several  attempts  have  been  made  to  accelerate  the  iterative 
convergence  of  this  method.  They  include  local  time  stepping,  implicit 
residual  smoothing,  enthalpy  damping  and  multigrid  techniques.  Also,  an 
extrapolation  procedure  based  on  the  power  method  and  the  Minimal  Residual 
Method  (MR1')  were  applied  to  the  finite  volume  Runge-Kutta  method.  In  the 
MRM,  a  weighted  combination  of  the  corrections  at  consecutive  iteration 
levels  is  extrapolated  and  the  weights  are  chosen  to  minimize  the  L2  norm 
of  the  future  residual.  The  extrapolation  was  performed  without 
considering  the  properties  of  the  governing  equations.  The  GNL.MR 
(Generalized  Non-Linear  Minimal  Residual)  method  utilizes  the  information 
about  the  governing  equations.  It  has  been  applied  successfully  to  a 
number  of  scalar  nonlinear  partial  differential  equations. 

Both  MRM  and  GNLNR  method  are  based  on  using  the  same  values  of 
optimal  weighting  factors  for  the  corrections  to  every  equation  in  a  system. 
Since  each  component  of  the  solution  vector  in  a  system  of  equations  has 
its  own  convergence  speed,  the  sequence  of  optimal  weights  could  be  allowed 
to  be  different  for  each  component.  This  concept  is  the  essence  of  the  DMR 
method.  Thus,  for  example,  we  combined  corrections  from  four  consecutive 
time  steps  by  introducing  four  weighting  factors  .0  each  of  the  four 
equations.  Hence,  a  set  of  sixteen  algebraic  equations  needs  to  be  solved 
to  determine  the  four  sequences  of  four  weighting  factors  in  each  of  them. 
The  DMR  method  requires  about  200%  more  storage  than  the  original 
non-accelerated  algorithm. 
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ABSTRACT 

A  new  algorithm  for  the  acceleration  of  explicit 
iterative  schemes  for  a  system  of  partial  differential 
equations  has  been  developed.  The  method  is  based  on 
the  idea  of  allowing  each  partial  differential 
equation  in  the  system  to  approach  the  converged 
solution  at  its  own  optimal  speed.  The  DMR 
(Distributed  Minimal  Residual)  method  allows  a 
separate  sequence  of  optimal  weighting  factors  to  be 
used  for  each  component  of  the  general  solution  vector. 
The  acceleration  scheme  was  applied  to  the  system  of 
time-dependent  Euler  equations  of  inviscid  gasdynamics 
in  conjunction  with  the  finite  volume  Runge-Kutta 
explicit  time-stepping  method  with  the  Jameson's 
Artificial  Dissipation  (AD)  terms  and  the  newly 
formulated  Physically  Based  Dissipation  (PBD)  model. 

The  PBD  model  uses  physical  dissipation  terms  from  the 
Navier-Stokes  equations  of  gasdynamics,  while 
enforcing  slip  boundary  conditions  of  inviscid 
gasdynamics  and  utilizing  spatially  varying  viscosity 
coefficients.  Tests  were  performed  for  various  flow 
conditions,  including  internal  flow,  flow  around  a 
cylinder  and  flow  over  an  airfoil  with  AD  and  PBD. 

Using  DMR,  between  30Z  and  70S  of  the  computational 
efforts  were  saved  in  the  subsonic  compressible  flow 
calculations . 


INTRODUCTION 

When  the  Euler  equations  of  inviscid  gasdynamics 
are  solved  using  a  central  difference  scheme  (e.g.,  a 
Runge-Kutta  time-stepping  scheme  [1]),  decoupling  of 
odd  and  even  grid  points  allows  oscillations  to 
develop  which  cause  instabilities  in  the  numerical 
algorithm.  These  oscillations  can  be  damped  by  either 
explicitly  or  implicitly  adding  a  certain  amount  of 
artificial  dissipation  (2). 

Contempc rary  artificial  dissipation  models  for 
central  difference  schemes  usually  consist  of  an  ad 
hoc  combination  of  second  order  and  fourth  order 
artificial  (non-physical)  dissipation  terms  [2].  The 
second  order  terms  are  used  to  damp  oscillations  in 
shock  regions,  while  the  fourth  order  terms  ensure 
"tonotonic  convergence  to  steady  state  in  smooth  flow 
regions  [1,3,4], 


Most  of  the  existing  artificial  dissipation 
formulations  are  intuitive  [5,1].  The  intuitive 
formulations  generate  artificial  dispersion  terms  (4j, 
which  are  partially  neutralized  by  adding  higher  order 
artificial  dissipation  terms.  Only  after  a 
trial-and-error  process  can  it  be  found  that  the 
coefficients  multiplying  second  and  fourth  order 
artificial  dissipation  which  are  appropriate  for  very 
low  speeds  are  orders  of  magnitude  smaller  then  the 
coefficients  that  are  appropriate  for  transonic  speeds. 
It  has  been  shown  that  using  different  amounts  of 
second  order  and  fourth  order  dissipation  can  produce 
different  numerical  results  that  are  often  misleading, 
especially  in  the  case  of  transonic  shocked  flows  with 
separation  [6,7] . 

It  can  be  concluded  that  the  intuitive 
formulations  for  artificial  dissipation  which  have 
been  favored  in  the  past  are  only  marginally  reliable. 
Their  accuracy  is  still  an  open  question  [6,7,8,9,10] 
since  there  is  no  known  exact  solution  to  the  Euler 
equations  for  a  shocked  flow  with  inviscid  separation. 
Thus,  the  existing  artificial  dissipation  models  are 
subject  to  constant  modifications  [4,11,12,13]  i.. 
order  to  meet  the  requirements  posed  by  different  flow 
speed  regimes. 

One  objective  of  this  paper  is  to  introduc 
physically  consistent  [14]  model  for  the  dissipat  :  r> 
to  be  used  in  the  numerical  solution  of  Euler  an 
Navier-Stokes  equations.  The  other  objective  is  to 
introduce  a  new  concept  for  convergence  acceleration. 

Several  att'mpts  have  been  made  in  the  past  to 
accelerate  the  iterative  convergence  of  the 
Runge-Kutta  method  [15].  They  include  local  time 
stepping  [1],  implicit  residual  smoothing  [1], 
enthalpy  damping  [1]  and  multigrid  techniques  [12,  Jo]. 
Also,  the  extrapolation  procedure  based  on  the  power 
method  and  the  Minimal  Residual  Method  (MRM)  were 
applied  [16]  to  the  Runge-Kutta  method.  In  the  MRM 
[16],  a  weighted  combination  of  the  corrections  at 
consecutive  iteration  levels  are  extrapolated  and  the 
weights  are  chosen  to  minimize  the  L2  norm  of  the 
future  residual.  Since  the  extrapolation  was 
performed  without  considering  the  properties  of  the 
governing  equations,  it  may  upset  the  solution 
procedure.  The  GNLMR  (Generalized  Non-Linear  Minimal 


Residual)  method  (17,18,19,20]  used  the  information 
from  the  governing  equations.  It  has  been  applied 
successfully  to  a  number  of  single  nonlinear  partial 
differential  equations  including  the  Euler  equations. 

Both  MRM  and  GNLMR  method  use  the  same  sequence 
of  optimal  weights  for  the  corrections  to  every 
equation  in  a  system.  Since  each  component  of  the 
solution  vector  has  its  own  convergence  speed,  the 
sequence  of  optimal  weights  should  be  allowed  to  vary 
from  component  to  component.  Thus,  the  objective  of 
this  paper  is  to  present  the  theory  constituting  the 
Distributed  Minimal  Residual  DMR  method  and 
todemonstrate  the  advantages  of  the  new  algorithm  with 
a  number  of  computational  examples. 

EULER  EQUATIONS  OF  CASDYNAM-CS 

The  two-dimensional  Euler  equations  in 
conservative  form  and  cartesian  coordinates  can  be 
written  as 


Q  ♦  E  ♦  F  -0 
t  x  y 

Here,  the  subscripts  t,  x,  y  represent  partial 
derivatives  with  respect  to  time,  and  to  x,  y 
coordinates,  respectively.  The  general  vectors 

Q,  E,  and  F  are  defined  as 


where  p,p,u,v,eQ  and  hQ  are  the  non-dimensional  values 
of  local  density,  thermodynamic  pressure,  x-component 
of  velocity,  y-component  of  velocity,  total  mass- 
specific  energy  and  total  mass-specific  enthalpy, 
respectively. 

Equation  of  state  for  a  calorically  perfect  gas 
can  be  expressed  as 


p  .  <Y-1)  (peo  -  I  ♦  ieiL))  (3, 

where  Y  represent*  the  ratio  of  specific  heats.  The 
total  mass-specific  enthalpy,  hQ,  is  defined  as 


For  the  analysis  of  flows  about  arbitrary 
geometries,  the  formulation  can  be  generalized  by 
u*ing»  s«y,  fixed  body-fitted  non-orthogonal 
coordinates  £  and  p,  so  that 


P  pu  pv 

Pu  PuU.ptJ  puV  .  pn 

q  ■  5  pv  e  -  -  pvu.pt  1  F  *  5  pw  +  Pny  (7) 


D  .  det  .  I  .  '  ^ 

a(x’y)  det  (§1^1  ly(  yn 


ex/D  *  V  VD  ‘  *  yt  :  V°  *  1y/D  -  (9) 

The  contravar iant  components  U,V  of  the  velocity 
vector  are  related  to  the  cartesian  components  u.v  as 
fol lows 


<x  S' 

%  V 


EXISTING  ARTIFICIAL  DISSIPATION  MODEL 

A  typical  stage  of  the  multistage  Runge-Kutta 
I  15]  time-stepping  scheme  for  the  Euler  equations 
11.4)  is 

..a'-1’.,  <„, 

where  Ptj  is  the  artificial  dissipation  [1,4 J  given 

as 

Pi,j'(Pl.l/2,j*  Pl-l/2,j)  *  ,Pi,j.l/2-  ri,j-l,2)  tl2) 

Calculation  of  the  artificially  dissipative  terms 
is  done  similarly  (1,4)  for  all  conservation  laws. 

For  example, 


P?  .  .«>  0t 

i+6/2,j  i+fi/2,j  ^i*5/2,j 


(4)  pf  cc 

'  Ci.6/2,j  (6i.6,j  °i.6,j  ~  °i , j  , j >  <I3) 

p  1  .  _(2)  nq 

i , j.ft/2  ri,j*6/2  Vi,j*4/2 

'  Ei  , j*4/2  (8i,j*6  Qi  J-4  ‘  6i,j  Qi^j)  (l4> 

where  4  -  ±1.  The  remaining  terms  are  defined  as 
foil ovs : 


$(x,y>  ;  R  -  q(x,y) 


Thus,  in  the  computational  (£,q) 
dimensional  Euler  equations  in  a 

domain,  the  tvo- 
strongly  conservative 

®i,j  -  — ; — 

(15) 

form  become 

Qi-4/2,j  ’  Qi*«,j  ‘  Qi,j 

(  16) 

Qt  4  V  Fn  *  0 

(6) 

9i+4,j  ‘  Qi*24,j  '  2Qi*S,j  *  Qi ,  j 

(17) 

where 

(18) 

QiJ  *  '  2Qi,j  *  Vl.j 
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with  similar  expression*  for  the  other  terms  of  this 
type.  The  coefficients  of  second  and  fourth  order 
artificial  dissipation  are  defined  (4,1], 
respectively,  as 


Ci+4/2,j  "  9i,j  * 


.) 

l.J  1  +  4. J 


Ci+4/2,j  ‘  n,ax  (°:  k’<4>'  ,t<2,max  ;  vi+6,j))  (20) 

with  similar  expressions  for  the  other  terms  of  this 
type,  where  K<2)  «  1/4  to  1/2  and  -  1/128  to  1/64 

are  typical  constants  |3).  Here,  the  local 
"directional  pressure  sensor"  is  defined  as 


ex  M  q 


•  Lnx  "yj  Lq; 


Here,  the  components  of  the  non-dimen6ional  viscous 
stress  tensor  expressed  ?n  terms  of  £,t)  coordinates 


pim  ~  2pi.j :  iizu 
pi+i.j  +  2pi.i  *  pt-i.j 


Similarly 


t  ■  p"  (C  u.  ♦  i)  u  )  ♦  K  (£  vP  ♦  i)  v  )  (27) 

xx  x  q  xi)  y  s  y  i) 

i  ■  p  ur  a  t)  u  ♦  £  ♦  n  v  )  (20) 

xy  y  C  y  n  x  £  'x  n 

and  the  non -dimension a  1  heat  conduction  flux  is 


^  lpi.J+1*2pi.i‘pi.j-1| 

PHYSICALLY  BASED  DISSIPATION  ( PBD)  MODEL 

Instead  of  using  an  intuitive  non-physical 
formulation  for  the  artificial  dissipation,  we  suggest 
that  the  dissipation  should  be  based  on  actual 
physical  dissipation,  that  is,  it  should  be  physically 
consistent.  We  propose  that  to  solve  the  Euler 
equations  of  inviscld  flow,  one  should  actually  solve 
the  complete  Navier-Stokes  equations  of  viscous  and 
heat  conducting  flow  subject  to  perfect  slip  boundary 
conditions  and  spatially  varying  coefficients  of 
viscosity  l 14].  Thus,  the  PBD  model  represents  a 
physically  consistent  formulation  since  the  Euler 
equations  of  inviscid  gasdynamics  represent  an  extreme 
case  of  Navier-Stokes  equations  when  the  physical 
dissipation  becomes  negligible. 

The  Navier-Stokes  equations  of  unsteady,  viscous, 
laminar  flow  allowing  for  heat  conduction  (assuming 
Fourier's  law)  expressed  in  non-dimensional  form  and 
non-orthogona 1  curvilinear  coordinates  can  be 
summarized  as 


q - K-__  (£  T  ♦  n  T  )  (29) 

(Y-1)M_  Pr  5  n 

Here,  y"  •  2y  ♦  K  is  the  longitudinal  viscosity 

coefficient,  is  the  Mach  number  of  the  uniform  flow 
at  infinity,  Pr  is  the  Prandtl  number  and  T  is  the 
absolute  temperature.  Since  Rankine-Hugoniot  shock 
jump  conditions  are  possible  only  (21]  if  Stokes 
hypothesis  (A/y  -  -2/3)  is  enforced,  we  use  this 
relation  in  actual  computations. 

In  the  PBD  formulation,  the  shear  viscosity 
coefficient,  y,  is  forced  to  vary  throughout  the 
flowfield  by  means  of  an  appropriate  "sensor".  The 
physical  thermodynamic  pressure,  p,  appears  in  the 
equations  of  gasdynamics  in  the  form  of  its  first 
derivative.  Consequently,  we  have  decided  to  use  the 
pressure  sensor  which  is  based  on  the  streamwise  first 
derivative  of  the  pressure,  that  is, 


>i.j  ’  (u2\-Vl72  (UPC  +  V  *  C  PsP 


Qt  *  Ee  *  Fn  -  S  <Ec  *  y 


where  R  is  the  Reynolds  number  and  Ee ,  F  incorporate 

e  t,  t) 

physically  dissipative  terms  due  to  shear  viscosity, 
secondary  viscosity  and  heat  conductivity.  The 
generalized  viscous  flux  vectors  are 


Here,  C  is  a  user  specified  constant .  Using 
numerical  experimentation,  we  have  found  that  10  <  C  < 

20  for  the  range  of  freestream  Mach  numbers  0. 1  <  M  < 

30.  We  have  experimented  with  a  number  of  different 
"sensors"  and  found  that  the  three-point  average 
streamwise  first  derivative  of  pressure  gives  a  robust 
scheme 

Vi  *  y+i.j +  "i.j  *  "i-i,j)  1 3  o,) 

Obviously,  this  is  just  one  among  many  possible 
suggestions  for  the  "sensor."  Other  choices  might  be, 
for  example, 

sensor  based  on  divergence: 


Vj  -  c I  7  • 


T  t 

yx  yy 


sensor  based  on  Mach  number: 


j-MCfMjO+M^V)/!  9| 


3  2- 
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sensor  based  on 


u  .  -  u«M^*M2V)/|  V| 
i.J  t  h 


DISTRIBUTED  MINIMAL  RESIDUAL  ( DMR )  METHOD 

Clobal  residual  of  the  finite  volume  method  at  time 
level  t  can  be  expressed  in  two  dimensions  as 

rt  -  fl  I? dS  -  il  (ff  *  ds 

where  S  is  the  surface  area  of  the  single  grid  cell 
and  components  Q,  E  and  F  of  the  generalized  solution 
vector  are  defined  in  Eq .  2.  In  the  case  of,  say, 
four-sten  explicit  Runge-Kutta  algorithm  one  needs 
four  intermediate  time  steps  to  advance  the  solution 
from  the  global  time  level  (t)  to  (t+1). 

We  plan  to  use  corrections  from  M  previous  consecutive 
time  levels  to  update  the  value  of  Q  to  (t+1)  global 
time  leve  1  .  Thus , 


Rt+1,  it  is  necessary  to  use  the  values  of 


Ql  *  I  n" 


“l  *1 
m  .  m 
U2  *2 
;  M  ;  M 
“L  6L 


and  are  the  corrections  for  each  of  the  t  *  1,...,L 

equations  in  the  system  (Eq.  2)  at  each  of  the 
m«l , . . . ,  M  global  time  steps.  Therefore,  substituting 
Eq .  36  in  Eq.  35,  the  new  local  residual  for  the 
single  cell  will  be 

rt  +  1  —  Jj  [|^  E  (Qc  *  2  *  |-  F  ( <5t  *  l  nm> )  dS  (18) 

m  m 

Using  a  Taylor  series  expansion  truncated  after  the 
first  term  results  in 


t+1  1  r  r r  .SE  m.  3  ,dF  nm.  ,  , . 

r  *  r  -  MI  l3?  (3Q  "  5  *  Si  %  "  )!  dS  (39> 

m 

Define  the  global  residual  as  a  sum  of  the  squares 
of  the  local  residuals,  that  is, 


Rl  *  1  l  (rC)*  (rC)  (AO) 
1  j 

where  I  and  J  define  the  grid  size  and  the  superscript 
*  ^.‘-’gnates  the  transpose.  Then,  the  global  residual 
at  the  next  global  time  level  will  be 


t  +  1  t*  r  /  t  r*  rr.3  . 3E  n .  3  ,3F  _n., 

R  *  l  l  (r  ~  £  lll3t  (3q  n  )f  a,,  (3q  n  >1  ds7 

l  j  n 

•  (rl  -  I  f|  ig  (g  nm)  *  Ij  -s) 


To  minimize  Rl  J 
«(  that  satisfy 


for  al)  m  and  f.  Thus,  from  Eq.  A1  and  Eq.  A2  it 
follows  that 
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where  - —  »  (A™  6,  )  and  6,  .  is  the  Kronecker  delta. 
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Notice  that 
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Note  that  A™  is  not  a  function  of  u's.  Then,  Eq .  (A3) 
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resulting  in  a  system  of  LxM  equations  for  the  unknown 

m 

distributed  optimal  acceleration  factors  .  In  tne 

case  of  two-dimes iona l  Euler  equations,  L  3  4.  Thus, 
we  must  solve  the  following  system  of  4xM  equations  in 

order  to  determine  the  4xM  optimal  values  of  . 
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We  have  decided  to  use  M  *  4,  that  is,  four 
consecutive  global  time  steps. 

Thus,  four  sequences  of  four  optimal  values  of  u 
were  used  in  Eq .  36  and  Eq.  37  to  update  the  solution 
to  the  next  global  time  level. 

RESULTS 

The  PBD  and  DMR  concepts  were  applied  to  three 
model  test  cases:  external  flow  around  a  cylinder, 
internal  channel  flow  past  a  101  circular  arc,  and 
external  flow  around  a  NACA  0012  airfoil. 

Figure  1  shows  fbe  65x33  0-type  compu t at i ona 1 
grid  around  a  cyiindet  .  1 igures  2  and  J  show  the 

convergence  histories  using  the  existing  Artificial 
Dissipation  (AD)  model  and  the  Physically  Based 
Dissipation  (PBD)  model  with  and  without  the 

application  of  DMR  with  *  0.2.  Using  DMR,  the 

number  of  iterations  n**eded  to  achieve  the  same  level 
of  residual  is  reduced  by  almost  5'1*.  The  savings  in 
computational  tine  is  about  "OX  (Figures  4  and  5) 
using  the  AD  and  the  PBD  model  with  the  application  of 

DMR  at  M  *  0.2.  The  savings  in  cpu  time  can  be  seen 

in  Table  l,  which  presents  the  run  times  and  residuals 
for  several  of  the  test  cases.  In  Table  1,  ResO  is 
the  starting  residual  and  R*»s  is  the  final  residual. 

Figures  6  and  7  show  convergence  histories  using 
the  AD  and  the  PBD  with  and  without  the  DMR 

model  for  the  M  =  0.4.  This  is  a  clear  indication 
that  the  present  formulation  of  DMR  is  incapable  of 


CPU  time  savings  as  the  locally  sonic  flow  conditions 
are  approached. 

At  M  =0.1  the  compressible  Euler  equations 
become  a  stiff  system  of  partial  differential 
equations.  Figure  8  shows  convergence  history  using 
the  AD  model  with  and  without  DMR.  The  large  CPU  time 
savings  demonstrate  the  ability  of  DMR  to  treat  stiff 
systems  of  equations  .(Figure  9). 

Figure  10  illustrates  the  65x17  H-type  chcmnel 
grid  with  a  101  circular  arc  bump  on  the  floor. 

Figures  11-12  show  that  at  ■  0.5  using  the  Au  -i - 
the  PBD  model  with  the  DMR  yields  about  5 OX  savings 
in  CPU  time.  Figures  13  and  14  show  that  using  AD  or 
the  PBD  model  with  DMR  at  M_  =  0.6  saves  less  than  50% 
in  CPU  time.  The  pressure  contours  using  the  AD  and 
the  PBD  models  were  identical  at  M^  *  0.6  (Figure  15). 

Figure  16  shows  the  65x33  C-type  clustered  grid 
around  a  NACA  0012  airfoil.  From  Figure  17  it 

appears  that  using  the  AD  model  and  DMR  at  M_  *  0.63 
does  nor  accelerate  convergence.  From  Figure  1ft,  the 

it  Is  clear  that  using  the  PBD  model  allows  the  DMR 
to  perform  better  even  for  this  transonic  shocked 
flow  case  resulting  in  over  30%  saving  in  the  CPU. 

Finally,  the  PBD  model  was  compared  to  the  AD 
model  by  applying  them  to  lifting  and  nonlifting 
transonic  flows.  The  1 2°x3 3  C-type  grid  around  a  NACA 
0012  airfoil  is  shown  in  Figure  19.  Figure  20  shows 

isobars  using  the  AD  model  at  =  0.8,  a  =  0.0°. 
Figure  21  shows  isobars  when  using  the  PBD  model  at 
M  *  0.8,  a  =  1.25°.  Again,  the  PBD  model  yields  a 
sharp  shock. 


CONCLUSIONS 

A  new  physically  based  dissipation  model  has  been 

presented.  Advantages  of  the  new  model  include: 

1.  The  second  order  dissipation  used  in  the  PBD 
model  represents  actual  physically 
consistent  dissipation  from  the  Navier- 
Stokes  equations  for  compressible,  viscous, 
heat  conducting  fluid  flow. 

2.  The  PBD  model  does  not  contaminate  the 
continuity  equation. 

3.  The  PBD  form*  at  ion  maintains  high  accuracy. 
Actually,  for  flows  with  stronger  shocks, 
the  PBD  formulation  gives  results  comparable 
to  TVD  schemas. 

4.  The  PBD  concept  can  be  applied  to  Navier- 
Stokes  equations,  too.  The  higher  order 
physically  consistent  dissipation  terms  can 
be  based  on  dissipation  due  to  radiation  heat 
transfer  and  heat  generation  due  to  chemical 
react  ions . 

5.  An  Euler  solver  with  the  PBD  formulation  easilv 
converts  to  a  Navi rr -Stokes  solver  by  fixing  ’he 
value  of  viscosit v  coefficient  and  by  specifving 
no-s 1  ip  boundary  conditions  or  solid  surfaces. 
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A  conceptually  new  method  termed  Distributed 
Minimal  Residual  ( DMR >  has  been  developed  and 
successfully  applied  to  the  acceleration  ot  an 
explicit  iterative  algorithm  for  the  numerical 
solution  of  a  nonlinear  system  of  Euler  equations 
governing  inviscid  gasdynamics.  The  main  idea  of 
using  a  separate  sequence  of  optimal  acceleration 
factors  for  each  of  the  equations  in  the  system  was 
theoretically  formulated  a  numerically  proven  on  a 
number  of  test  cases.  This  means  that  the  partial 
differential  equations  governing  mass,  x-momentum, 
y-momentum  and  energy  conservation  were  accelerated 
according  to  their  own  separate  optimal  sequences  of 
acceleration  factors  that  have  a  common  objective  of 
minimizing  the  global  residual  of  the  entire  system  at 
each  consecutive  integration  time  step.  DMR  in  its 
present  form  works  best  for  low  Mach  number  flows  when 
the  Euler  equations  become  exceedingly  stiff. 
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Table  i.  Circle;  Computational  O-type  grid  (6$  x  33) 

Comparison  of  CPU  time  (sec)  for  Artificial 
Dissipation  (AD)  and  Physically  Based 
Dissipation  (pbd) 
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Circle:  Convergence  history;  AO;  1%,-  0.2. 
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Circle:  Convergence  history;  PBD;  M„-  0.2. 


87 


NACA0012:  Convergence  history;  PBD;  M®»  0,53 


NACA0012:  Computational  grid  !129  x  33) 


Figure  20. 


Figure  21 
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ABSTRACT 

A  new  algorithm  for  the  acceleration  of 
iterative  schemes  for  the  numerical  solution  of 
systems  of  partial  differential  equations  has  been 
developed.  The  method  is  based  on  the  idea  of 
allowing  each  partial  differential  equation  in  the 
system  to  approach  the  converged  solution  at  its 
own  optimal  speed.  The  DMR  (Distributed 
Minimal  Residual)  method  allows  a  separate 
sequence  of  optimal  weighting  factors  to  be  used 
for  each  equation  in  the  system.  The  acceleration 
scheme  was  applied  to  the  system  of  time-dependent 
Euler  equations  of  inviscid  gasdynamies  in 
conjunction  with  the  finite  volume  Rational 
Runge-Kutta  (RRK)  explicit  time-stepping  algorithm. 
Using  DMR  without  mul tigri ddi ng,  between  30i  and 
70S  of  the  total  computational  efforts  were  saved 
in  the  subsonic  compressible  flow  calculations. 

DMR  method  in  its  present  version  seems  to  be 
especially  suitable  for  stiff  systems  of 
equations.  It  required  less  than  double  amount  of 
storage  of  the  original  non-accelerated  algorithm. 


INTRODUCTION 

One  of  the  successful,  explicit  methods  usee 
to  solve  Euler  and  Navi er-Stokes  equations 
governing  compressible  flows  subject  to  the 
various  flow  conditions  is  the  Rational 
Runge-Kutta  (RRK)  time-stepping  algorithm  [1,2]. 

It  is  based  on  the  finite  volume  technique  with 
2nd-**th  order  blended  non-physical  (artificial) 
dissipation  [1].  Several  attempts  have  been  made 
to  accelerate  the  Iterative  convergence  of  this 
method.  They  include  local  time  stepping  [1], 
Implicit  residual  smoothing  [1],  enthalpy  damping 
[1]  and  multigrid  techniques  [33.  Also,  an 
extrapolation  procedure  based  on  the  power  method 
and  the  Minimal  Residual  Method  ( MRM )  were  applied 
[3]  to  the  finite  volume  Runge-Kutta  method 
together  with  multigridding.  In  the  MRM  [3],  a 
weighted  combination  of  the  corrections  at 
consecutive  iteration  levels  is  extrapolated  and 
the  weights  are  chosen  to  minimize  the  Lj  norm  of 
the  future  residual.  The  extrapolation  was 
performed  without  considering  the  specific 
properties  of  the  governing  equations.  The  GNLMR 
(Generalized  Non-Linear  Minimal  Residual)  method 
[*i,5,6, 7]  utilizes  the  information  from  the 


governing  equations.  It  has  been  applied 
successfully  to  a  number  of  scalar  linear  and 
nonlinear  partial  differential  equations. 

Both  MRM  and  GNLMR  method  use  the  same  values 
of  optimal  weights  for  the  corrections  to  every 
equation  in  a  system.  Nevertheless,  since  each 
component  of  the  solution  vector  in  a  system  of 
equations  has  its  own  convergence  speed,  the 
sequence  of  optimal  weights  could  be  allowed  to 
vary  from  equation  to  equation.  The  authors 
believe  that  this  concept  underlying  the 
Distributed  Minimal  Residual  (DMR)  method  is 
similar  to  the  general  idea  behind  the 
preconditioning  techniques.  With  the 
preconditioning,  the  eigenvalues  of  the  system  are 
changed  so  that  the  different  CFL  (Courant- 
Friedrichs-Levy)  number  can  be  used  for  each 
characteristic  variable.  This  paper  presents  the 
theory  constituting  the  DMR  method  and 
demonstrates  the  advantages  of  the  new  algorithm 
with  a  number  of  computational  examples. 
Applications  of  the  DMR  to  the  system  of  Euler 
equations  of  inviscid  gasdynamies  are  presented. 
The  formulation  can  be  equally  well  applied  to 
other  systems  of  differential  equations  and  to 
otner  types  of  numerical  integration  algorithms. 


TIME-DEPENDENT  EULER  EQUATIONS  OF  INVISCID 
GAS  DYNAMICS 

The  system  of  time-dependent  Euler  equations  of 
gasdynamies  in  two-dimensional  space  can  be 
written  in  a  general  conservative  form  as 


3Q  ^  3E  +  3F 
3t  36  3n 


(1) 


where  the  global  solution  vectors  combining  mass, 
x-mcmentum,  y-roomentura  and  energy  conservation 
equations  are  defined  as 
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Tne  terms  on  tne  rignt  hand  siaes  of  c.q.  &  are 
similar  [1].  For  example, 
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Here,  p,  u,  v,  p,  e0  are  the  density,  x  and  y 
components  of  the  velocity  vector,  thermodynamic 
pressure,  and  mass-specific  total  energy, 
respectively.  In  addition,  U,  V,  E,  n  ana  p  are 
the  contravari ant  velocity  vector  components, 
non-orthogonal  curvilinear  computational 
coordinates,  and  determinant  of  the  Jacobian 
transformation  3(E ,  n)/3(x,y ) ,  respectively . 

The  contravarlant  components  U  and  V  of  the 
velocity  vector  in  the  body-conforming  (E.n) 
coordinate  system  are  given  by 
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where  the  second  and  fourth  order  coefficients 
multiplying  the  flux  Derivative  terms  are  flow 
adaptive  coefficients.  Tne  scaling  with  the  area 
b  anti  the  local  time  step,  it,  is  included  [S]  to 
correspond  to  the  formulation  of  the  Euler 
equations  In  tne  transformed  plane.  A  pressure 
sensor  is  introduced  to  locate  regions  requiring 

large  amounts  of  artificial  dlsrlpatlor..  It  is 
based  on  the  second  derivative  of  pressure  [1,9] 


U  -  E  u  ♦  E  v 
x  y 
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The  total  energy  per  unit  mass  for  a  calorically 
perf ect  gas  Is 
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The  flow  adaptive  coefficients  are  then  calculated 
[1]  as: 


-  V  +  (5) 

where  cv  Is  the  specific  heat  at  constant  volume 
:  and  T  Is  the  absolute  temperature.  The 
determinant  of  the  Jacobian  geometric 
i  transformation  matrix  is 
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FINITE  VOLUME  RUNGE-KUTTA  TIME-STEPPING  ALGORITHM 

In  the  finite  volume  method  [ 1 D ,  the 
governing  equations  are  Integrated  over  each 
computational  cell  In  the  (E.n)  computational 
plane.  Kith  the  help  of  the  divergence  theorem, 
the  surface  Integral  is  transformed  into  a  sum  of 
line  Integrals.  These  integrals  are  discretized 
with  the  assumption  that  the  fluxes  are  constant 
along  the  cell  faces.  Each  quantity  at  the  cell 
face  is  evaluated  as  the  average  of  the  values  at 
the  neighboring  cell  centers  (cell  centered 
scheme) . 

The  cell  centered  finite  volume  method  is 
Identical  to  the  central  difference  scheme  on  a 
uniform  grid.  It  Is  known  that  the  central 
difference  scheme  produces  odd-even  decoupling. 

To  suppress  this  tendency, the  artificial 
d.ssipation  terms  are  added  to  the  discretized 
equation  [1].  The  mixture  of  2nd  and  Ath  order 
artificial  dissipation  terms  [1]  was  used. 


where  d  is  the  artificial  dissipation  operator  and 
Q  is  the  vector  defined  In  Eq.  2.  The  two  terms 
on  the  right  hand  side  of  Eq.  7  are  contributions 
from  the  two  computational  directions.  They  can 
be  written  [1 ]  as: 


The  system  of  time-dependent  Euler  equations 
is  known  to  be  of  hyperbolic  type  and  the  boundary 
conditions  should  be  applied  according  to  the 
direction  of  the  characteristics.  At  the  inflow 
and  outflow  boundaries,  the  incoming  Riemann 
Invariant  is  specified  and  the  outgoing  Riemann 
Invariant  Is  extrapolated  from  the  Interior  points. 
Also,  the  entropy  and  the  tangential  velocity  are 
prescribed  at  the  inflow.  At  the  outflow,  these 
quantities  are  extrapolated  from  the  interior  of 
the  domain. 

At  the  solid  wall,  the  normal  momentum  equation  is 
used  to  evaluate  the  wall  pressure.  The 
contravarlant  velocity  component  U  at  the  ghost 
cells  inside  the  solid  body  Is  extrapolated,  while 
the  oontravariant  velocity  component  V  is 
reflected  from  the  wall. 

An  explicit  Runge-Kutta  time-stepping  [2,1] 
scheme  Is  used  to  evolve  the  solution  in  time. 

The  Ath  order  Runge-Kutta  scheme  is  given  by 


1-1/2J’ 


dlj-M/2  *  dij-1/2 


- 

Qn 

(12) 

c(1)  - 

Qn  - 

at 

A 

C  N  C C  0  ? 

-  d5(0)) 

03) 

■ 

CVI 

wc / 

Qn  - 

at 

3 

(NC0) 

-  dC(0)) 

(U) 

Q^)  • 

Qn  - 

At 

2 

(NQ(2) 

-  dQW) 

03) 

- 

c  - 

At 

-  dC(0>> 

06) 

q"*1  - 

07) 

2 


where  N  is  the  discretization  operator  or  tne 
finite  volume  method.  The  artificial  dissipation 
is  evaluated  at  the  beginning  of  each  Runge-Kutta 
global  step  and  it  was  not  updated  during  the 
Intermediate  steps.  Linear  stability  analysis 
indicates  that  the  explicit  Runge-Kutta  scheme  is 
stable  if  CFL  £  2.S.  If  the  grid  spacing  in  (C.n) 
plane  is  uniform  if  *  in  «  i,  then  the  time  step 
is  given  [9]  by 


Using  a  Taylo--  series  expansion  truncated  after 
the  first  term  results  in 
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Define  the  global  residual  fit  as  a  sum.  of  the 
squares  of  th<  local  residuals,  that  is, 


where  a  is  the  local  speed  of  sound  and  the 
combined  time  step  [S]  is 
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where  I  and  J  define  the  grid  size  and  the 
superscript  *  designates  the  transpose  of  an  array, 
Then,  the  global  residual  at  the  next  global  time 

level  will  be 
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residual  of  the  finite  volume  method  at 
time  level  t  can  be  expressed  as 
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where  S  is  the  surface  of  the  single  grid  cell  and 
components  0,  E  and  F  of  the  generalized  solution 
vector  are  defined  in  Eq.  2. 

We  plan  to  use  corrections  from  M  consecutive  time 
levels  to  update  the  value  of  Q  to  (t-1  )  global 
time  level .  Thus , 
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T c  minimize  R**' ,  it  is  necessary  to  use  the 
values  of  uT  that  satisfy 


for  all  m  and  1.  Thus,  from  Eq .  29  and  Eq .  26  it 
follows  that 
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and  Ll  are  the  corrections  and  are  the  weights 

for  each  of  the  ...,l  equations  in  the  system 
(r,q.  2)  at  each  of  the  m  consecutive 

globa*  time  levels.  Therefore,  upon  substituting 
Eq.  22  in  Eq.  21,  the  new  local  residual  for  the 
single  cell  will  be 
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Notice  that 
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and  4  t  is  the  Kronecker  delta. 
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Note  that  A^  is  not  a  function  of  u's.  Then,  Eq . 
( 31  )  becomes 
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resulting  in  a  system  of  LxK  equations  for  the  LxK 

unknown  optimal  acceleration  factors  u.®.  In  the 
case  of  two-dimensional  Euler  equations,  L  *  A. 
Thus,  we  must  solve  simultaneously  the  following 
system  of  AxK  equations  in  order  to  determine  the 
AxK  optimal  values  of 
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We  have  decided  to  combine  four  consecutive 
tine  step3  (M-A).  Since  the  twodimensional  Euler 
equations  form  a  system  that  has  four  equations  (L 
-  a),  these  four  sequences  of  four  optimal  values 
of  u  can  then  be  used  in  Eq.  23  and  Eq.  22  to 
update  the  solution  to  the  next  global  time  level 
t*1 . 


The  matrices  9E/90  and  3F/9C  that  are  needed  for 
evaluation  of  the  coefficients  in  the  above  matrix 
are  given  as: 
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whe-e  h  is  the  specific  enthalpy  per  unit  rr.ass  and 
Y  in  the  ratio  of  specific  heats  for  a  caioricaliy 
perfect  gas. 


In  addition  to  the  computer  memory  required 
by  the  original  nor.-accele-ated  scheme  [1], 
additional  memory  is  needed  to  implement  the  DMR . 
If  the  grid  points  are  IxJ  and  we  use  M  global 
consecutive  time  levels  to  update  the  solution, 
then  for  the  twc-dlmer.sional  problem  the  extra 
memory  requirement  is  approximately  L  x  (2-*-M)  x 
(1-2)  x  (J-2)  and  for  the  three-dimensional  Euler 
equations  the  extra  memory  requirement  is 
approximately  L  x  (3*M)  x  (1-2)  x  (J-2)  x  (K-2). 
In  the  two-dimer.si onal  case  thir.  represents 
approximately  1 505  increase  and  in  the 
three-dimensional  case  this  represents 
approximately  175'  increase  in  memory  requirement 
over  the  original  nor-accelerated  [1]  algorithm. 

Three  different  methods  were  tested  for  the 
boundary  conditions  on  the  residuals  in  the 
integrals  of  Eq.  35.  The  first  method  was  to  set 
the  residuals  at  the  ghost  cells  to  be  zero.  The 
second  method  calculates  the  residuals  at  the 
ghost  cells  from  the  boundary  conditions.  The 
third  method  extrapolates  the  residuals  from  the 
Interior  of  the  flowfield.  It  was  found  that  the 
third  method  gives  the  best  results. 


All  computations  were  performed  on  a  VAX 
l  ‘  /'6550  computer  In  a  single  precision  mode.  Tr.e 
f:rst  secjer.ce  cf  tests  cf  DidR  was  performed  or. 
the  internal  two-dimensional  (L-4)  flow  problems 
by  combir  'g  four  consecutive  global  time  steps 
(M-4).  This  means  that  a  1  6x  1  6  matrix  (  Eq  .  4",) 
needs  to  be  inverted.  Figure  1  snows  the 
computational  grid  for  a  10'  thick  circular  bump 
1  r.  a  two  dimensional  channel.  Tne  grid  size  is 
c'cxl?  points.  The  calculations  were  started  with 
uni  form  flow  and  the  DMR  was  applied  once  afte- 
every  30  Iterations.  Figure  1  shows  the 
convergence  histories  of  subsonic  flow 
calculations  with  K*,-  0.5.  The  number  of 
iterations  needed  to  achieve  the  same  level  of 
residual  is  reduced  almost  t  60'.  The 
convergence  with  the  DMR  shows  smaller 
oscillations  than  that  of  tne  original  [l]  scheme. 
It  is  expected  that  this  behavior  continues  to  the 
machine  accuracy.  The  saving  in  computational 
time  Is  about  50%  for  this  test  case. 

The  constant  pressure  contour  plots  of  the 
entire  flow  field  for  both  non-aceelerated  finite 
volume  RP.E  scheme  and  DMR  accelerated  finite 
volume  RRK  scheme  are  shown  in  Fig.  4  and  Fig.  5, 
respectively.  The  difference  between  the  two  is 
not  discernable  lr.  these  contour  plots  thus 
confirming  that  DMR  method  does  not  adversely 
influence  the  quality  of  the  solution. 

Results  of  the  second  test  case  are  presented 
in  Figs.  6,7,8  and  9.  The  entire  flow  field  is 
subsonic  with  M„  *  0.55.  For  this  test  case,  the 
saving  was  almost  40%  In  CPU  time.  It  is 
noticeable  that  the  convergence  history  shows  more 
oscillatory  behavior  than  for  the  case  with  M  - 
C.5.  Another  subsonic  (Men*  0.6)  test  case  was 
tested  and  the  results  are  shown  in  Figs.  10,  11, 
12  and  13  demonstrating  that  a  considerable  amount 
of  computation  effort  was  saved. 

Figs.  14  and  15  show  the  convergence 
histories  for  the  transonic  shocked  flow  case  with 
Mcd  ■  0.675  which  is  less  than  the  flow  choking 
Mach  number  of  this  channel.  Results  indicate 
that  with  the  DMR,  the  convergence  rate  is  not 
improved. 

Similar  trends  were  observed  when  solving 
Euler  equations  for  a  flow  around  a  circle.  An 
0-type  grid  consisting  of  64x32  grid  cells  was 
used.  For  a  moderately  compressible  subsonic  flow 
(Mas-  0.3),  DMR  saves  (Figs.  18  and  19) 
approximately  45%  of  CPU  time.  It  generates 
results  (Fig.  20)  that  are  practically 
indistinguishable  from  the  non-accelerated  scheme. 
When  the  critical  free  stream  Mach  number  M<t>  -  C.4 
was  used.  Fig.  21  indicates  and  Fig.  22  confirms 
that  the  DMR  method  In  its  present  form  offers 
practically  no  gain  when  compared  with  the 
non-accelerated  algorithm  although  the  computed 
surface  Mach  numbers  (Fig.  23)  are  equally 
accurate.  Thus,  both  Ni's  bump  cas*.  and  circle 
case  indicate  that  DM?,  method  in  its  present 
f c—uiati on  offers  no  advantages  at  transonic 
speeds.  On  the  other  hand,  the  system  of  Euler 
equations  becomes  stiff  as  the  Mach  number 
decreases,  thus  rapidly  reducing  the  convergence 
rate  of  tne  non-accelerated  scheme.  When  using  M* 

=  C.l  (an  almost  incompressible  flow)..  Figs.  24 
ar.d  25  demonstrate  that  DM?,  offers  over  70% 
savings  in  the  CPU  time  over  the  non-accelerated 
scheme.  Fig.  26  indicates  difference  in  the 
computed  su-face  Cp  values  after  1200  iterations. 
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In  orde-  do  account  for  the  different  local 
characteristic  behavior  of  the  transonic  flow,  it 
should  be  possible  to  use  different  sets  o. 
weights  for  different  regions  of  the  flow-field. 
Also,  t.ne  artificial  dissipation  terms  could  be 
incorporated  in  the  formulation  of  the  SMR.  Ir. 
addition,  the  optimal  frequency  of  applying  the 

DMP  needs  t-e  investigated.  In  the  present 
investigatic  ,  DMR  was  applied  by  combining  four 
consecutive  time  steps  after  every  thirty  time 
steps . 

Notice  that  all  numerical  results  were 
obtained  without  tne  standard  acceleration 
techniques  such  as  explicit  and  implicit  residual 
smoothing,  enthalpy  damping,  multigriddlng  and 
vert  or i zat i on .  These  methods  could  be  added  to 
further  accelerate  the  algorithm. 

CONCLUSIONS 

A  conceptually  new  method  termed  Distributed 
Minimal  Residual  (DMR)  has  been  developed  and 
successfully  applied  to  the  acceleration  of  an 
explicit  finite  volume  iterative  algorithm  for  the 
numerical  solution  of  a  nonlinear  system  of  Euler 
equations  governing  lnviscid  gasdynamlcs.  The 
main  idea  of  using  a  separate  sequence  of  optimal 
acceleration  factors  for  each  of  the  equations  in 
the  system  was  theoretically  formulated  a 
numerically  demonstrated  with  a  number  of  examples. 
This  means  that  the  partial  differential  equations 
governing  mass,  x-momentum,  y-momentum  and  energy 
conservation  were  accelerated  according  to  their 
own  separate  optimal  sequences  of  acceleration 
factors  that  have  a  common  objective  of  minimizing 
the  global  residual  of  the  entire  system  at  each 
global  time  level.  The  method  seems  to  offer 
significant  time  savings  especially  for  stiff 
systems  cf  differential  equations. 
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Figure  1.  Computational  grid  for  a  10'  thick 

circular  arc  airfoil  on  the  bottom  wall 
of  a  straight  two-dimensional  channel 
(Ni's  bump). 
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Figure  2.  Comparison  of  convergence  rates  in 
terms  of  Iteration  numbers:  non¬ 
accelerated  and  DMR  accelerated 
algorithm  for  Ki's  bump  with  t'a,  -  0.5. 


"Ul/i 


Figure  if.  Constant  pressure  contours  f or-  non¬ 
accelerated  algor: the,  for  Ill's  Pump 
with  Men'  0.5. 
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Figure  3.  Comparison  of  convergence  rates  in 

terms  of  the  CPU  time:  non-accelerated 
and  DMR  accelerated  algorithm  for  Nl's 
burp  with  Mo,-  0.5. 
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igure  5.  Constant  pressure  contours  for  DMR 
accelerated  algorithm  for  Nl's  bump 
with  K«,.  C.5. 
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Figure  6.  Comparison  of  convergence  rates  in 
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accelerated  and  DMR  accelerated 
algorithm  for  Mi's  bump  with  Keo  *  0.55. 
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INTRODUCTION 


{ 


One  of  the  successful,  explicit  methods  used  to  solve  Euler  and  Navier-Stokes 
equations  governing  compressible  flows  is  the  finite  volume  Runge-Kutt3  time-stepping 
algorithm  [1].  Several  attempts  have  been  made  to  accelerate  the  iterative  convergence 
of  this  method.  These  acceleration  methods  are  based  on  local  time  stepping  [1],  implicit 
residual  smoothing  [1],  enthalpy  damping  [1]  and  multigrid  technique;  [2).  Also,  an 
extrapolation  procedure  based  on  the  power  method  and  the  Minimal  Residual  Meihcn. 
(MRM)  were  applied  f2]  to  the  Jameson’s  multigrid  algorithm.  The  MRM  has  not  beet 
shown  to  accelerate  the  scheme  without  multigridding.  It  uses  same  values  of  optimal 
weights  for  the  corrections  to  every  equation  in  a  system,  i  each  component  of  the 
solution  vector  in  a  system  of  equations  is  allowed  to  have  its  own  convergence  speed, 


then  a  separate  sequence  of  optimal  weights  could  be  assigned  to  each  equation.  This 
idea  is  the  essence  of  the  Distributed  Minimal  Residual  (DMR)  method  [3]  which  is  based 
on  the  General  Nonlinear  Minimal  Residual  (GNLMR)  concept  [4], 

TIME- DEPENDENT  EULER  EQUATIONS 

The  system  of  time-dependent  Euler  equations  of  gasdynamics  in  a  two-dimensional  space 


can  be  wrirten  in  a  genera!  conservative  form  [1]  as 

§  +  f  +  f  =°  ( 

rt  or? 

where  the  global  solution  vectors  combining  mass,  {-momentum.  r?-momentum  and 


energy  conservation  equations  are  defined  as 


P  U 

puU  + 
h'U  + 
+ 


P  3’ 

ArV  +  r;xp 
pvV  +  T)  p 

Ae0  +  p)V 


Here,  t,  p,  u,  v,  p,  ec  are  time,  density,  x  and  y  components  of  the  velocity 
vector,  thermodynamic  pressure,  and  mass-specific  total  energy,  respectively.  In  addition. 
U,  V,  r,  and  D  are  the  contravariant  velocity  vector  components,  non-orthogonal 
curvilinear  computational  coordinates,  and  determinant  of  the  Jacobian  of  the 
transformation,  3(f„v)  '(^x,y),  respectively. 

DISTRIBUTED  MINIMAL  RESIDUAL  (DMR)  METHOD 

Local  residual  of  the  fini'  -volume  method  at  the  global  time  level  t  can  be  expressed  a? 

=  //  f¥ -  -if  (f +  % >  dS  (3, 

where  S  is  the  surface  of  the  single  grid  cell  and  the  components  Q(,  El  and  Fl 

are  defined  in  Eq.  2.  In  the  DMR,  corrections  from  M  consecutive  time  levels  are  used 

to  update  the  value  of  Q  to  (t+1)  global  time  level.  Thus, 

Qt+1  =  Q*  +  I  (4) 

m 

where 


Here,  Arr‘  are  the  corrections  computed  with  the  basic  algorithm  and  u™ 

are  the  weights  for  each  of  the  1=1 _ L  equations  in  the  system  (Eq.  2)  at  each  of  the 

m=l,...,  M  consecutive  time  levels.  Therefore,  the  new-  local  residual  for  the  single  grid 
cell  will  be 


(2) 


Define  the  global  residual  R1*1  at  the  global  time  level  (t+1)  as  a  sum  of  the 
squares  of  t!  •  local  residuals,  that  is. 


Rt+1  =  V  V  (rt  +  1)"  (rt41) 


vhere  1  and  J  define  the  grid  size  and  the  superscript  *  designates  the  transpose. 


objective  is  to  find  optimum  v.lues  of  L  sequences  of  M  values  of  up  t 
minimize  the  global  residual  Rt+1  at  the  next  global  time  le  el  (t+1).  To 
it  is  necessary  to  use  the  values  that  satisfy 

=  o 

'rr 

(Tj 

1 

for  all  m  and  1.  Thus,  from  Eq.  7,  8  and  9  it  follows  that 


that  will 


minimize 


(r*  -  J  fj  [|  <g  W  -  |  tf)]  dS>’ 
.  i  rr  ,1  ^  +  *  (£El  dSi  =  o 


fja-  =  (A%) 


and  6k,  is  the  Kronecher  delta.  Notice  that 


c*  _  y  s£L  un 
V  q 


i  JJ  l*  (aQt  aup>  a?  a,p)J 


Note  that  Am  is  not  a  function  of  u/s.  Then,  Eq.  (10)  becomes 

i 


(13) 


I  J  M  L  ,  1 

YYYY  ( An )  A"^  =  E  £  ( r!)  Am  (U> 

ryYY  «  *  q  y  j  > 

Let 

i  J  .  i  J 

cnm  =  v  V  ( Ar  )  Am  Bn‘  =  v  v  ( r*  j  Am  (15) 

ql  T  y  *•  1  !  i  J  1 

The  result  is  a  system  of  LxM  equations 

y  cnm  +  tn  Cnm  +  un  cnm  + . +  ^  cnm)  =  Bm  (16) 

n  il!  :  21  3  31  L  LI  I 

for  the  LxM  unknown  optimal  acceleration  factors  v™.  The  DMR  applied  to  the 

finite  volume  scheme  [1]  in  two-dimensional  case  needs  approximately  150%  increase  and 
in  the  three-dimensional  case  it  needs  approximately  175%  increase  in  computer  memory 
over  the  original  non-accelerated  algorithm  [1].  Boundary  conditions  on  the  residuals  in 
the  integrals  of  Eq.  13  used  extrapolation  of  the  residuals  from  the  interior  of  the 
flowfield. 

RESULTS 

All  computations  were  performed  on  a  VAX  11/8550  computer  in  a  single  precision 
mode.  The  first  sequence  of  tests  was  performed  on  the  interna)  two-dimensional  (L  =  T 
flow  problems  by  combining  four  consecutive  time  steps  (M  =  4).  This  means  that  a 
16x16  matrix  (Eq.  16)  needs  to  be  inserted.  The  test  geometry  was  a  10%  thick  circular 
half  airfoil  on  a  wall  of  a  straight  two  dimensional  channel.  The  H-npe  grid  size  was 
65x17  points.  The  calculations  were  started  with  uniform  flow  3nd  the  DMR  was  applied 
once  after  every  30  steps  performed  with  the  original  unaccelerated  code  [1].  Figures  1 
and  2  depict  the  convergence  histories  of  flow  calculations  with  Mw  =  0.5  and 
=  0.675.  For  the  entirely  subsonic  flow  ( M'30  =  0.5)  the  number  of  iterations  needed  to 
achieve  the  same  level  of  residual  is  reduced  almost  by  60%,  while  the  saving  in 


T 


computational  time  is  about  5QV  Both  figures  indicate  that  DMR  in  its  present  version 
does  not  accelerate  transonic  flow  =  0.675)  computations.  Superimposed  constant 
pressure  contours  (Fig.  3)  of  the  entire  flow  field  for  both  the  non-accelerated  and  the 
DMR  accelerated  schemes  confirm  that  DMR  method  does  not  adversely  influence  th 
quality  of  the  solution. 

The  secc.  d  test  case  was  a  flow  ar  und  a  circle.  An  O-tvpe  radially  c<  stered  grit, 
consisting  of  64x32  grid  cells  was  used.  We  applied  DMR  after  every  60  iterations  by 
combining  four  consecutive  time  levels.  When  the  critical  free  stream  (M^.  =  0.4)  was 
used,  Figs.  4  and  5  indicate  that  the  DMR  method  in  its  present  form  offers  practically 
no  gain.  At  very  low  free  stream  Mach  numbers  the  system  of  Euler  equations  becomes 
very  stiff,  thus  rapidly  reducing  the  convergence  rate  of  the  non-accelerated  scheme. 

On  the  other  hand,  when  using  Ma  =  0.1,  the  DMR  offers  over  70%  savings  in  the 
CPU  time  (Fig.  5)  over  the  non-accelerated  scheme. 

Notice  that  all  numerical  results  were  obtained  without  the  standard  acceleration 
techniques  such  as  explicit  and  implicit  residual  smoothing,  enthalpy  damping, 
multigridding  and  veetorization.  These  methods  could  be  added  to  further  accelerate  the 
algorithm.  The  method  seems  to  offer  substantial  time  savings  when  applied  to 

compressible  flow  codes  at  low  Mach  numbers. 

CONCLUSIONS 

A  new  method  for  the  acceleration  of  explicit  iterative  algorithms  for  the  numerical 
solution  of  systems  of  partial  difierential  equations  has  been  developed.  The  method  is 
based  on  the  idea  of  allowing  each  partial  differential  equation  in  the  system  to  approach 
the  converged  solution  at  its  own  optimal  speed  while  at  th'  same  time  communicating 
with  the  rest  of  the  equations  in  the  system.  The  D MR  (Distributed  Minimal  Residual) 
method  computes  a  separate  sequence  of  optima!  acceleration  factors  to  be  used  for  each 


r? 


component  of  the  general  solution  vector.  The  acceleration  scheme  was  applied  to  the 
system  of  tinie-dependent  Euler  equations  of  inviscid  gasdynamics  in  conjunction  with  the 
finite  volume  Runge-Kutta  explicit  time-stepping  algorithm.  Using  DMR  without 
multigridding.  between  30%  and  %>%  of  the  total  computational  efforts  were  saved  in  the 
subsonic  compressible  flow  calculations.  The  DMR  method  seems  to  be  especially  suitable 
for  stiff  systems  of  equations  and  can  be  applied  to  other  systems  of  differential 
equations  and  other  numerical  algorithms. 
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Figure  1.  Convergence  histories  in  terms  of  iteration  numbers:  non-aecelerated 

( - land  DMR  accelerated  ( -  )  algorithm:  channel  flow.-. 

Figure  2.  Convergence  histories  in  terms  of  the  CPU  time:  non-accelerated  ( - i 

and  DMR  accelerated  algorithm:  channel  flow. 

Figure  3.  Constant  pressure  contours  for  DMR  accelerated  algorithm:  transonic 

channel  flow  with  MWi  =  0.675. 

Figure  4.  Convergence  histories  in  terms  of  iteration  numbers:  non-accelerated 

( - )  and  DMR  accelerated  ( _ 1  algorithm:  circle  flow 

Figure  5.  Convergence  histories  in  terms  of  the  CPU  time:  non-accelcrnted  ( - ) 

and  DMR  accelerated  ( _ )  algorithm:  circle  flow 


6  / 


VTM VT3TTt>  O.T  Tm^^r 
i\ui— — »r\  v-  -  —  u.  *-»-  j.  — 


<£p- 


